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We study excitonic effects in two-dimensional massless Dirac fermions with Coulomb interactions 
by solving the ladder approximation to the Bethe-Salpeter equation. It is found that the general 
4-leg vertex has a power law behavior with the exponent going from real to complex as the coupling 
constant is increased. This change of behavior is manifested in the antisymmetric response, which 
displays power law behavior at small wavevectors reminiscent of a critical state, and a change in 
this power law from real to complex that is accompanied by poles in the response function for finite 
size systems, suggesting a phase transition for strong enough interactions. The density-density 
response is also calculated, for which no critical behavior is found. We demonstrate that exciton 
correlations enhance the cusp in the irreducible polarizability at 2/cf, leading to a strong increase 
in the amplitude of Friedel oscillations around a charged impurity. 

PACS numbers: 71.45.Gm,73.22.Pr,71.10.-w 



I. INTRODUCTION 

Graphene is a two-dimensional honeycomb lattice of 
carbon atoms. Near zero doping, the low energy quasi- 
particle states are well-described by two-dimensional 
massless Dirac fermions (MDF's). Graphene supports 
two species of these, centered at the two inequivalent 
corners of the Brillouin zone. At low or zero doping, the 
properties of graphene are in many ways quite different 
than those of a doped semiconductor, in large part be- 
cause of the unusual properties of MDF'si^. 

One interesting class of questions about graphene 
involve Coulomb interactions. In a seminal study, 
Gonzalez, Guinea and Vozmediano2 demonstrated that 
weak Coulomb interactions in undoped graphene are 
marginally irrelevant, invalidating the basic premise of 
strong screening that underlies the standard treatment 
of an electron g weakly interacting Fermi liquid. 

Moreover, simple estimates of the strength of Coulomb 
interactions, characterized by an effective fine structure 
constant j3 — e 2 /ehvF, where vp is the speed of electrons 
near a Dirac point and e is the effective dielectric constant 
due to a substrate upon which the graphene may be ad- 
sorbed, suggest that Coulomb interactions are effectively 
large {(3 > 1) if e ~ 1. This suggests that properties 
of graphene near zero doping should be rather different 
than those of non-interacting MDF's; for example, a gap 
may open in the quasiparticle spectrum^, so that rather 
than behaving as a metal the system would be insulat- 
ing. In real experiments there is little evidence for such 
dramatic effects of interactions, most likely because dis- 
order effects overwhelm those of Coulomb interactions 2 . 
If so, one may suppose that interaction effects could be- 
come apparent if sufficiently clean graphene samples can 
be created. In this work, we study this theoretical clean 
limit, and focus on unusual properties which can emerge 
in linear response functions for two-dimensional MDF's 
due to Coulomb interactions. As we shall see, these have 



important consequences for the induced charge distribu- 
tion around an impurity, and suggest that the MDF de- 
scription breaks down even before a gap opens in the 
spectrum. 

The unusual effects of a Coulomb potential for MDF's 
are already apparent in their response to a Coulomb 
charge Ze, even when there are no interactions among the 
electrons themselves. The wavefunctions for this problem 
are essentially exactly calculable 5 - - — . One finds that the 
m-th circular component of an electron wavefunction has 
the short distance form ip m (r) ~ rV^+Va) 2 -^ 2 ^ at 
a distance r from the impurity. This is an unusual situ- 
ation in that the exponent is a function of the impurity 
charge, so that the wavefunctions have a non-analytic 
dependence on the potential strength at short distances. 
This cannot be reproduced at any finite order in per- 
turbation theory. When Z/3 is above a critical value, the 
short distance exponent becomes complex, corresponding 
to a maximal penetration of the centrifugal barrier (the 
effective potential due to the angular momentum in the 
radial equation for the wavefunction) by the electrons. 
We refer to this phenomenon as "Coulomb implosion" , 
and it is analogous 5 to the breakdown of the vacuum in 
the vicinity of a highly charged nucleus in QED— . In 
graphene, this breakdown is accompanied by the forma- 
tion of a charge cloud around the impurity with density 
falling off as 1/r 2 , which is absent for Z/3 below its crit- 
ical value. The propagation of a short distance effect 
(penetration of the centrifugal barrier) to long distances 
(appearance of the charge cloud) is one of the special 
properties of MDF's, and it reflects the absence of any 
length scale in the Dirac equation itself. As we shall see, 
there are many-body analogs of these phenomenon which 
become apparent in some of the linear response functions. 

To address the many-body problem, it is preferable to 
assess the non-analytic content of a linear response func- 
tion as a function of momentum rather than position. To 
see how this might be done, we revisit the problem of non- 
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interacting MDF's in the presence of a Coulomb impurity, 
and analyze how the short distance behavior described 
above is manifested in a momentum representation. Not 
surprisingly, we find that the scattering wavefunctions, 
when expressed in a momentum representation, display 
power law behavior at large momentum, with an expo- 
nent that changes from real to complex for Zf3 exceeding 
the same critical value found in the real space analy- 
sis. The momentum space analysis in terms of scatter- 
ing states naturally suggests that one might find simi- 
lar behavior in vertex functions evaluated in the ladder 
approximation^. Ladder diagrams play an important 
role in interacting systems because they allow one to in- 
corporate excitonic correlations in virtual particle-hole 
pairs that are generated by the interaction. We analyze 
this approximation for a generic four-point vertex func- 
tion of MDF's with Coulomb interactions, and find both 
the non-analytic power law behavior at large momentum 
and a change from real to complex exponent, in this case 
when (3 exceeds some critical value. The analysis suggests 
the system undergoes a quantum phase transition at this 
critical value, since the exponent necessarily behaves in 
a non-analytic way as a function of the parameter /3. 
Interestingly, we shall see that the equations for the ver- 
tex function suggest that this transition is infinite order, 
suggesting the transition may be in the same universality 
class as classical two dimensional systems undergoing a 
Kosterlitz-Thouless (KT) transition^. 

To relate the four point vertex function to measurable 
quantities, we use it to form three-point vertex func- 
tions which can then be used to directly compute lin- 
ear response functions. Two are of particular interest. 
The density response function (equivalently, the density- 
density correlation function) expresses the screening re- 
sponse of the system to an external potential, including 
due to impurities or inhomogeneities in the system. We 
find that the non-analyticity of the four-leg vertex does 
not present itself in this particular quantity. Neverthe- 
less, we demonstrate that exciton effects, as expressed in 
the ladder diagrams, have important quantitative effects 
for for doped systems, where we find a strong enhance- 
ment of Friedel oscillations around a charged impurity 
relative to the non-interacting case. 

The other important response function involves the 
charge imbalance between the two sublattices, the an- 
tisymmetric response function. (Some results on this 
were reported by us previously^.) Here we indeed find 
non-analytic behavior analogous to the Coulomb impu- 
rity problem for non-interacting MDF's. Although this 
non-analyticity is apparent at large wavevectors in the 
four-leg vertex, the absence of length scale in the problem 
leads to power law behavior emerging at small wavevec- 
tors in the response function. In particular this applies 
to an impurity placed asymmetrically with respect to the 
graphene sublattices, which is known to induce very dif- 
ferent charge responses on them 15 - 18 . The sublattice- 
antisymmetric component of this acquires a power law 
tail, with exponent that changes from real to complex 



above a critical value of /?. Interestingly, when a finite 
size cutoff is included in the calculation, we find poles in 
the response function, suggesting a quantum phase tran- 
sition to a state with different charge densities on the 
sublattices, and hence a state with spontaneously bro- 
ken chiral symmetr y 19 ' 20 . These poles however merge 
together into a branch cut in the thermodynamic limits, 
suggesting a state fluctuating among ones with the chi- 
ral symmetry broken in different possible ways, and no 
mean-field mass gap. This may represent a phase transi- 
tion that is a precursor to one in which a real gap develops 
in the spectrum^. 

This article is organized as follows. We begin with a 
study of the single impurity problem in the momentum 
representation in Section |TT] and identify the signatures 
for the short-distance power law behavior and Coulomb 
implosion. With this simpler example to guide us, in Sec- 
tion [TTT] we study the ladder approximation to the Bethe- 
Salpeter equation for the 4-leg vertex. This approxima- 
tion is the many-body analog of the Lippmann-Schwinger 
equation studied in Section fll] We show that the ver- 
tex function has a power law behavior in momentum, 
and the exponent changes from real to complex above 
a certain value of coupling constant /3. The interesting 
properties of the antisymmetric response are discussed 
in Section HVl Finally, in Section [V] we focus on the 
density-density response function in the ladder approxi- 
mation, and demonstrate the important quantitative ef- 
fects caused by excitonic correlations. 



II. IMPURITY PROBLEM IN THE 
MOMENTUM REPRESENTATION 

We begin by discussing the problem of MDF's in the 
presence of a Coulomb impurity in terms of scattering 
states. The standard (Lippmann-Schwinger) equation for 
scattering states 2 ^ takes the form 

^+\x)=^°\x) + / d 2 yG DE (x~y)V(y)^+Hy), (1) 



where G DE is the (matrix) Green's function deter- 
mined by the differential equation (DE) (—hvpa ■ p + 
el)G UE {x) = SW(x), V(y) = Ze 2 /e\y\ is the impurity 
potential and ip( ' is an eigenstate for V = 0. Fourier 
transforming Eq. (TTJ, we have 

(2) 

It is convenient to express the Lippmann-Schwinger equa- 
tion in terms of angular momentum states. Introducing 
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the angular components 
Jo 27r 



(3) 



P\) = J] e- in ^- e ^f n (p/p')/p', (4) 



(5) 



where Op is the angle between the vector p and the x axis, 
and using 



G DE (p) = 



1 



e 2 — (hv F p 2 ) 



[el + Hv F ( Px a x + Py a y )}, (6) 



where a x and a v are the Pauli matrices, Eq. ([2j may be 
written in the form 



4 + i = v£(p) 



(7) 



Ze 2 p 



(hv F p) 2 



+hv F p I xf^ {m+ i ) {x)^ n+1 {xp)dx 



Ze 2 p 



(0) 



(8) 



[fepp / x/_ m (a;)'^ 1 t^(icp)drc 



e 2 - (hv F p)- j 

/•QQ 

+e / x/_ (m+ i)(x)^;^ +1 (xp)dx]. 
Jo 

Motivated by the observation that power law behav- 
ior emerges in the wavefunctions at small distances, we 
search for power law solutions at large wavevector. Using 
the ansatz 



^a + l(p) 



P 



for p — >• oo , 



(9) 



and neglecting terms of lower order of p, Eqs. [7] may be 
written as 



1 

Zf3I m (s) 



-Zj3I m+ i(s) \ ( a 



where 



Im(s) = X S f- m (x)dx. 

Jo 



(10) 



(11) 



Eqs. [10] will have non- vanishing solutions provided 

l-(Z/3) 2 I m (s)I m+1 (s) = Q. (12) 

One may easily show that I m (s) has a minimum at 
s = 3/2 for any integer to, so for Zf3 larger than a crit- 
ical {Zf3) c , the solution s to Eq. [T2] becomes complex. 
For m = 0, (Z0) c = 1/2; for m = 1, {Z0) c = 3/2, etc. 
This change of behavior corresponds to that found in the 
real space analysis of the Coulomb impurity problem-,. 
For Zj3 > [Z0) c (for a given to), the complex values of s 



Cause tpa^m (p) to oscillate at small r, with no well-defined 
value of ipatrn(r) as r —t 0. This leads to an ill-defined 
problem unless a boundary condition for small but fi- 
nite r is imposed 5 . This suggests that some quantities 
are sensitive to the short scale cutoff in the problem, a 
behavior which we will see holds true as well for some 
response response functions when Coulomb interactions 
are included. Using the above analysis as a guide, we 
now turn to this more complicated problem. 



III. 



GENERAL 4-LEG VERTEX 



As we saw in the last section, the signature of Coulomb 
implosion in the momentum representation is that the 
exponent of the power law of the wavefunction becomes 
complex. The basic physics in Eq. (Q} is clear: in a 
perturbative expansion in V, the electron can be scat- 
tered arbitrarily many times by the impurity, and the 
non-analytic, power law behavior emerges from a su- 
perposition of all these possibilities. If we substitute 
the electron-impurity scattering with electron-hole scat- 
tering due to Coulomb interactions, we may expect a 
many-body analog of both the power law behavior and 
of Coulomb implosion. The signature of these should 
be contained in the general 4-leg vertex function in the 
electron-hole channel. 

To compute this, we notei^ that the ladder approxima- 
tion to the 4-leg vertex (Fig. [JJ has the same structure of 
multiple scattering of an electron from a hole as does an 
expansion of Eq. [TJ in powers of V. The Bethe-Salpcter 
equation resulting from this ladder sum has the form 12 



r a p,-/6(Pl,P2]P3,P4,) = U(Pl -P3)5ajSps 



(13) 



where T is the vertex function whose arguments are three 
momenta [spatial components (momentum) p and time 
component (frequency) po ], U(p) — 2~Ke 2 /t\p\ is the 
Coulomb interaction, and 



G io) (po,P) 



(14) 



Po + A 4 /^ + v fP ■ <? 



{po + pi/h) 2 - (v F p) 2 + iS ■ sgn(p )(po + fi/K) 

is the time-ordered Green's function for non-interacting 
Dirac fermions, and we have allowed the possibility of a 
non-zero chemical potential /i. 

Following Ref. I12L we introduce two new functions Q 
and x y i a the relations 



T a p n s(pi,P2;P3,P4) 

d 3 q 



(15) 



U(q)Q a p,ys{pi - q,P2 - q;P3,P4), 



(2tt) 3 

Xafi n s{p,p';P) = (16) 
/ ^Qa^siP + \p,p- ~P;p' + \p,p' \p)- 
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FIG. 1: Ladder approximation to the Bethe-Salpeter equa- 
tion. 



In these expressions, P represents the four-momentum of 
the particle-hole pair, which may be understood as en- 
tering the vertex at the bottom of the diagrams in Fig. 
[T] and exiting at the top. In this interpretation, p then 
represents the relative momentum of the pair before the 
collision, and p' the momentum afterwards. These equa- 
tions allow Eq. (fl~3|) to take the form 



d 2 q 



(17) 



+K a p. f _ l „(p, P) 



(2ir} 



where 



K a p,^(p,P) 



d ^l {p+ \p)G[%- l -P). 



,( 18 ) 

Note that because we have integrated over po, there is no 
frequency dependence in K or x^- The quantity x can 
be used to directly to construct static response functions, 
as we shall see below. 

In the case of the Dirac particle colliding with a 
Coulomb impurity we found power law behavior for large 
momenta, specifically from collisions involving a large 
change in momentum (\p\ » [p'l in Eq. [5J) Since in 
the two-body problem, the particle and hole scatter from 
one another, we search for analogous behavior in the ver- 
tex functions at large momentum difference \p — p' \ . For 
example, for fixed P and p' , when \p\ — > oo, the equa- 
tions for xn,7<5 an d X22,75 become 



Xu,ys(p) = 

1 f d 2 q 



(19) 



4hv F \p\ J (2tt) 2 
1 f d 2 q 



U{p- q)\xn^s(q) - X22, 7 s(q)} 



(20) 



Ahv F \p \ J (2tt) 2 



U(p- q)[-Xu,-ys(q) +X22,~,s(q)} 



Note in these expressions terms of order P/p, p' jp have 
been dropped, and p and P are suppressed in the argu- 
ments of x- 

Introducing circular moments as before, i.e. 



r^e-^x(pl, (21) 



and using the ansatz 



(m) , v _ Wyyj 
Xii n s\P> - pS ■ 

(m) 
22,75 
pS 



X2I %(P) 



c. 



(22) 
(23) 



the m-th angular component of the above coupled equa- 
tions reduces to 



p -^m(^) -^m(s) \ ( C^l 



f = 0- (24) 



Nontrivial solutions to this set of homogeneous linear 
equations may be found if 



2//3 = I m (s), 



(25) 



which determines the exponent s for a given /3. The 
critical /? where the exponent of the power law becomes 
complex is (3 C — 2// m (3/2). We thus see the many-body 
problem has a Coulomb implosion instability analogous 
to what is found in the Coulomb impurity problem for 
non-interacting MDF's. A natural interpretation for this 
instability is that it indicates a transition to a gapped, 
excitonic insulator stat o 19 i 20 . However, as we shall see 
below, when analyzed in terms of the appropriate linear 
response function, such an interpretation is only consis- 
tent when considering a finite size system 1 ^; in the ther- 
modynamic limit the transition is likely to a state with 
a fluctuating mass gap. 

Xii\s an d X2T75 are not tne om Y components of the 
vertex function which display instabilities. An analogous 
instability appears for Xi2,7<5 a nd X21.75, which are gov- 
erned by the equations 



Xl2,j6(f) = 

1 f d 2 q 



4hv F \p\ J (2tt) 2 

X21, 1 s(f) = 

1 f d 2 q 



4hv F \p\ J (2tt) 2 



(26) 

U(p- q)[xi2,,s(q) - e^^i^l], 

(27) 

U(p - q)[~e 2ie ? X i2ns(q) + X2i, 7 «(5)]- 



The equations for the coefficients in the power law ansatz 
are 



■f — Im{ s ) Im+2{ s ) 



C 



m) 



Us) l-I m+2 (s))[ C ^p J" ' (28) 



so that the equation for s is 



4//3 = I m (s) + I m+2 {s). 



(29) 



The critical /3 for the 12 and 21 components of \ is 
p' c = 4/[J m (3/2) + J ro+2 (3/2)] > /3 C . This suggests that 
Xulyg an( ^ ^22*75 can develop complex exponents before 
Xi™ 7 5 an d X2i \s as ft IS increased from small values. 



5 



We will see more generally that certain combinations 
of the vertex functions do not appear to develop power 
law behavior at all; most importantly this appears to be 
the case for the vertex function relevant to the density- 
density response function. In this context, we note that 
the coefficients satisfy the conditions Cj^' s = —G'i" 1 ^ 



and cQls 



" C 2™7« > so that Xll,7«(P) 



y 22,7(5 

and Xi2,js{p) = — X2i,"fs{p)- These relations among the 
components of Xa0,-<s can a l so be seen from Eqs. [TTland 
IT51 when expressed in terms of circular moments. The 
density-density response turns out to involve Xaa,~n ( re ~ 
peated indices here are summed), so that the power law 
behavior is canceled away. 

We have verified our analysis by numerically solving 
the Bethe-Salpeter equation in the form of Eq. (TTT1) . We 
use polar coordinates for the integration and change each 
dimension of the integration to a discrete sum, indepen- 
dent of the other dimension, i.e. 



d0 



dqf(0,q) 



No 

i=i 



E 



WjWuqj), (30) 



where A is the momentum cutoff and / denotes a gen- 
eral integrand, and the discretization in each dimension 
is done according to the Gauss-Legendre rule 2 ^. Now the 
integral equation is changed to a set of linear equations 
and can be solved using any existing subroutine, e.g. the 
appropriate routine in the Lapack package. Some exam- 
ples of our results are shown in Fig. [2J As can be seen, 
the vertex functions follow power law forms, and the ex- 
ponents are close to those expected from our asymptotic 
analysis, both below and above the critical value of j3. 

An interesting aspect of these results is that the change 
of the exponent s = s' + is" from real to complex at a 
finite value of /3 (i.e., s" is zero for /3 < j3 c but is non- 
vanishing for /3 > p c ) suggests that quantities calculated 
from it will have a non-analyticity at /3 C . However, any 
such singularity must be of infinite order. For exam- 
ple, if there is a cusp in x at /3 = /3 C , then there would 
be a contribution of the form B a p n s(p,p'; P)d(/3 — j3 c ) 
in d 2 x/df3 2 . Differentiating Eq. (|I7[) with respect to /3 
twice, and requiring that the coefficients of 5(/3 — (3 C ) on 
both sides of the equation are the same, we find an equa- 
tion for B, 



d 2 q 
(2tt) 2 



(31) 

[U(q)\p=p c ]B flUt7S (p- q,p';P). 



This is nothing but a homogeneous version of Eq. (flTl) . 
with f3 [in U (q)] set to /3 C . If Eq. (|3Tj) has a non-vanishing 
solution, then Eq. (|17j) would also have a contribution 
from such a solution, and we would expect \ to be di- 
vergent at /3 = (3 C . Our explicit solutions, both in the 
asymptotic and the numerical analysis, show that this is 
not the case, so that no cusp can be present. Similarly, 
higher order derivatives of x a ls° cannot have cusps. It 



0.100 
0.050 

0.020 
0.010 
0.005 

0.002 
0.001 





FIG. 2: (Color online) Some examples of power laws from 
solving Eq. (|f 7[) numerically. P = and p' — fcp/2. (a) 
f3 — 0.1 < f} c , the solid lines are fits with the model C/p s , 
where p = p/kF- For SRxiiii we g e * s — 2.02 from fitting, 
s = 1.95 from solving Eq. ({2SJ. For KX1212 we S et s = 3-09 
from fitting, s = 3.05 from solving Eq. (b) /3 = 1 > j3 c . 

The solid line is a fit with C cos (s" log p + S) /p s . The fitting 
gives s' = 1.57 and s" = 0.78, while Eq. ^ gives s' = 1.5, 
s" = 0.60. 



follows that the singular behavior at /3 = /3 C - and any 
phase transition it may represent - is of infinite order, a 
property it shares with KT transitions^. In our analysis 
of the antisymmetric response below we shall sec further 
hints of a connection to the KT universality class in this 
system. 

The symmetries of the components of x discussed 
above led to a cancellation of the power law behavior in 
the density-density response function. This same sym- 
metry suggests that in a response function of the form 
a apXai3.-yS, with (T z the Pauli matrix, the power law 
should be retained. This represents the response of the 
density difference between the A and B sublattices due 
to a potential that is antisymmetric in sublattice index. 
Any potential that breaks sublattice symmetry will have 
a component of this antisymmetric response, so that it is 
in principle physically accessible. As we shall show next, 
the antisymmetric response does capture the power law 
as well as the Coulomb implosion physics. 
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FIG. 3: (Color online) (a) Diagrammatical equation for the 
3-leg vertex T^(k, q) with a z as the zeroth order vertex (the 
shaded cross in the figure), (b) Diagram for M(q). 



FIG. 4: (Color online) (a) An insertion for the reducible 
diagram for M(q). (b) The diagram equation for the 3-leg 
Coulomb vertex r a2( 9 2 (k, q). 



IV. ANTISYMMETRIC RESPONSE 



We define the antisymmetric response as 

M{q) = ~j dt([m z (-q,t),m z (q,Q)]), (32) 



where A is the area of the sample, repeated indices are 
summed, and henceforward we will set ft = 1. In princi- 
ple M(q) can be determined experimentally by measur- 
ing the screening charge induced by an impurity placed 
asymmetrically with respect to the two sublattices, i.e. 
anywhere except in the center of a hexagon, or at the 
middle point of a carbon-carbon bond. The difference 
between densities on the two sublattices rh z is an interest- 
ing operator because when < rh z (q) > is non- vanishing, 
there is a dynamically generated Dirac mas o 19 i 20 . 



A. Diagrammatic Expansion and Ladder 
Approximation 

The diagrammatic representation for M(q) is illus- 
trated in Fig. 3(b) along with the summation of ladder 
diagrams, Fig. 3(a)| representing our approximation for 
the vertex function. Notice there are no bubble diagrams 
in the diagrammatic expansion of M(q); there are only 
irreducible diagrams. Reducible diagrams for this quan- 
tity turn out to vanish, as we now show. Any reducible 
diagram will have at its end an insertion of the form il- 
lustrated in Fig. 4(a) in which there is a Coulomb vertex 
r. This insertion represents a multiplicative contribution 
to the diagram of the form 



d 3 k 
(2tt) 3 
d 2 k 



(2tt 
d 2 k 



<A G^ 2ai (k + q)Gf l02 (k)T a202 (k, q) (34) 

2 ^Qift-^aiftc^fe (&! 9)r Q2/ 3 2 (k, q) 



In these expressions the overhead tilde denotes quantities 
pertaining to 3-leg vertex. Note also that we have set 



qo = in the second line. Assuming the Coulomb vertex 
is given by a ladder sum, the equation determining x is 
[Fig. [4(b)] 



d 2 q' 



(35) 



UiWDKasPwnih ~ ^q)T lll2 {k - q,q), 



(27T) 2 

so that x satisfies 

Xaf}(k, q) = K afja2a2 (k, q) 
d 2 q' 



(36) 



-K a p a2 p 2 (k,q) 



(2nf 



U(W\)Xa 2 p 2 (k - cf,q). 



More explicitly, the 11 and 22 components of Eq. 
in the undoped case (n = 0), are 



Xn(k,q) 



d 2 k' 



1 



-l{O k -0+) 



+ 



1 



2(k+\k + q\) 2(k + \k + q\) 



(37) 



(27T) 



X22(k,q) = 



U(k-k') [xil(k',q)-e 
I _ e i(e k -e+) 



_ p -K8k-6+)~ 



X22(k',q) 



2(k+\k + q]) 2(k+\k + q}) 



(38) 



d 2 q' 



U(k-k') \- e ^ k -o + )^ l{ k' iq)+ X22 (k',q) 



where 9k and 6 + are the angular coordinates of the two- 
dimensional momenta k and k + q respectively. Without 
loss of generality, we can set q y = 0, from which it is easy 
to see that X22(k x , — k y ; q) = Xn(fc x , k y ; q). When this 
is substituted into the last of Eqs. [M] one readily sees 
that this insertion vanishes. Thus our approximation for 
M(q) includes only the irreducible ladder diagrams. 

The calculation of the M(q) follows steps very anal- 
ogous to those described for the 4-leg vertex, and were 
outlined in Ref. The equation for the 3-leg antisym- 
metric vertex [Fig. 3(a)| is 



f^(k,q)=a* a p 



-^U(W\)K a0lS (k-<f,q) 



T™ s (k-<f,q). 



(39) 
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Defining X^{k,q) = K a0lS (k,q)f^ s (k,q), one finds 

Xap( k ,q) = K af}l s{k,q)a z lS + K a ^ s (k,q) 
d 2 q' 



U(W\) X %(k-<f,q). (40) 



(2nf 



This quantity is related to the susceptibility by 

d 2 k 



M{q) 



(27T) 



(41) 



B. Solutions at Long Wavelengths 

In what follows we focus on the long wavelength limit 
(small q), so we drop all terms of 0(q 2 ) and higher. Using 
a circular moment expansion one finds 



fe'*'/o(x)4 (0) ( fc '^ ( 42 ) 



Here we used the superscript (0) to denote the circular 
component m = 0, and the underlined indices are not 
summed over. Note that we have introduced both an 
ultraviolet cutoff (A ~ 27r/a, a = lattice spacing) and an 
infrared cutoff (fco ~ 2n/L, L = linear size of system). 

Defining X M ^(k,q) = VppXpa (k,q), in the limit 
q — s- the solution to Eq. (j42|) may be written in the 
form x M (°)(fc, 0) = (a) > where F obeys the inte- 

gral equation 



pij -i 



(43) 



Note that F depends on the ratio fc/A, a reflection of the 
fact that the original Hamiltonian has no intrinsic length 
scale, so (in the limit ko — > 0) k can enter only in this 
ratio. For fc/A <C 1, one easily confirms that Eq. (|4"3")l is 
solved by a power law F (j) ~ (A/fc) s , with s going from 
real to complex above some critical j3. This is precisely 
the behavior we identified in the 4-leg vertex; unlike what 
one finds in the density-density response case, the power 
law is not canceled upon forming the 3-leg vertex from 
the 4-leg vertex. 

Eq. (|4"3")l may be readily solved numerically. For small 
/3, the solution is indeed a power law, provided k 3> fco 
[see Fig. 5(a) inset]. For large enough j3, the solution is 
consistent with a power law of complex exponent, such 
that F becomes oscillatory with a power law envelope 
[Fig. [5(a)]. Interestingly, M(q^ 0) = / jB?x(k,g^ 0) 
also has a series of divergences [Fig. 5(b)| . Formally, one 
may understand the occurrence of these poles by thinking 
of the solution in terms of the inverse of 1 — j3L, where L 
is the integral operator on the right hand side of Eq. (|43| . 
Divergences then occur as 4 crosses successive eigenval- 
ues of L. The presence of such poles suggests a phase 
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FIG. 5: Solutions of Eq. (gHJ) with k /A = 1CT 1U . (a) j3 = 0.5. 
Because the plotted is \F\, the oscillations appear as cusps. 
Note the amplitude of the oscillation scales roughly as 1/ \fk. 
Inset: F for /3 = 0.3. It is clearly a power law except for k 
close to fco ■ (b) The antisymmetric response M as a function 
of the interaction strength /3. 



transition into a state with a spontaneously generated 
M(q — > 0) becoming a Dirac mass, i.e., chiral symmetry 
breaking. However, the positions and weights of these 
poles are sensitive to fco , the infrared cutoff due to the fi- 
nite system size, and merge together in the L — > oo limit 
to introduce a branch cut in F as a function of (3. We 
discuss the significance of this below. 

For small but nonzero q, it is interesting to compute 
the correction AM (q) — M(q) — M(0). The equation for 
the corresponding AF has a form very similar to Eq. (|43l) , 
with only the "1" replaced by an inhomogeneous term, 
which is proportional to q 2 /k 2 for k q. The AM(q) 
resulting from this then vanishes with an exponent that 
varies with j3. The inset of Fig. [6] illustrates a typical re- 
sult for P not too large; the exponent as a function of /3 is 
illustrated in the main panel of Fig.[B] One physical con- 
sequence of this is that the difference in charge between 
sublattices for an impurity placed asymmetrically with 
respect to the sublattices will fall off with a /3-dependent 
power law at large distances, behavior which may be ob- 
servable with a local scanning probe. 

The result illustrated in Figs. [5] and [5] have a num- 
ber of interesting consequences for interacting electrons 
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FIG. 6: The exponent in AM(q) as a funciton of ft. Inset: 
AM(q) for ft = 0.3. 



in undoped graphene. For ft < ft c , we see that there 
are indeed power law correlations at long distances in 
quantities that are in principle measurable, with an ex- 
ponent varying continuously with ft. This means that the 
weak-coupling many-body groundstate possesses a basic 
property of a critical phase. For ft > ft c the exponent 
becomes complex, as in the noninteracting Coulomb im- 
plosion problem. In the interacting many-body case, the 
susceptibility M(q) of Eq. (|55t diverges for fc > 0. This 
strongly suggests a quantum phase transition to broken 
symmetry state with staggered charge orde r 19 ! 20 . 

The evolution of this system, from a state with power 
law correlations to one in which a transition occurs when 
this power reaches some limiting value, is highly remi- 
niscent of the phenomenology of systems which undergo 
a KT transition^. In such systems the transition in- 
dicates the appearance of a correlation length, which 
equivalently indicates that a gap develops in the exci- 
tation spectrum. This behavior in general should be sig- 
naled by a divergence in an appropriate response func- 
tion. However, the presence of many such divergences as 
a function of ft suggests there are different ways to break 
the symmetry. For a system with finite system size, we 
expect this chiral symmetry-breaking to occur as ft in- 
creases from small values in a way that is consistent with 
the first such pole. As we show below, the separation 
between neighboring poles vanishes only logarithmically 
as fco ~~ ► 0, so that finite size may in fact be important 
for realistic system sizes. 

Nevertheless, the merging of these poles suggests 
that something else must happen in the thermodynamic 
limits. The merging of these poles as fco results in a 
a continuous function with a branch point at ft c . We in- 
terpret this latter non-analytic behavior as the signal of a 
phase transition. Since it is a result of the merging poles, 
a natural interpretation is that the instability is into a 
state involving fluctuations among different realizations 
of a chiral order parameter which, if quiescent, would 
produce a gapped exciton phas o 19 i 20 . We speculate that 



with further increase in ft, one of these orderings could 
be favored over the others, resulting in a true condensed 
phase. This would be consistent with results of quantum 
Monte Carlo calculations^. 



C. An Analytical Model 

A fuller understanding of Eq. (I43[) may be arrived at 
with a model kernel of the form 



f (x) = 9(1 - x) + -0{x - 1). 

x 



(44) 



This has the same behavior as the real kernel at large and 
small x, and is simple enough to allow analytic solutions. 
We have verified numerically that the results for F and 
M are qualitatively very similar to those obtained with 
the correct / . 



1. Solution by Wiener-Hopf Method 

With this model kernel, assuming that the integral con- 
verges as q — > 0, which will be checked after the solution 
is found, we obtain, in the fco — > limit, the integral 
equation 



Let us change to more convenient variables via 



k 



-t K _ -t' 



A =e ' A =e 



(46) 



and rename the function for which we are solving, as well 
as the kernel, 

F (j)=9(t), fo(jj =K(t-t% (47) 
to obtain the rewritten integral equation 

oo 

9(t) -~j dt'e^'Kit - t')g{t') = 1. (48) 
o 

Note that physically meaningful values of t are non- 
negative. An important point is that if the limits of inte- 
gration over t' had been (— co, oo) one could have solved 
the equation trivially by Fourier transformation. Since it 
is over the half-line one has to use the more sophisticated 
Wiener-Hopf method 24 . 

One first extends the definition of g so that it has the 
entire real line for its domain, defining 



g(t) = g+(t) + g-(t), 



(49) 
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where g+(t) is nonzero only when t > and g~(t) is 
nonzero only for t < 0. As part of the solution one ob- 
tains both g + and g_. Defining R{t) = e t K{t) we can 
now extend the range of integration over (— oo, oo) as 
long as one integrates only g + : 



2 



dt' R(t-t')g + (t') = l-g-(t) 



r+(t)+r_(t) 
(50) 

where r+(t) = Q(t) and r_(t) = [1 - g-(t)] 9(-t). 

One can now solve the equation by Fourier transforma- 
tion. The crucial point is that since g+ is nonzero only for 
nonnegative values, if it vanishes as t — > oo, its Fourier 
transform has poles only in the lower half-plane of com- 
plex to, while <?_ has poles only in the upper half-plane. 
This gives us the extra information needed to solve for 
both g + and <?_. In general, there is no need for g + to 
vanish as t — > oo, which would correspond to the original 
function F vanishing as k — > 0. In fact, one expects F to 
diverge with a power law as k — > 0. To incorporate this 
expectation, we define g+{t) = e st h+(i), with h + vanish- 
ing as t — ¥ oo. The new function h + satisfies an integral 



equation with a modified kernel R s (t) 



M*)-| / dt'R 8 (t-t')h+(t') = e- st {r + (t) + r-(t)). 

— OO 

We now take the Fourier transform of both sides, using 
the explicit form of /o, which corresponds to 



R s (t) = e^-^Oi-t) + e- st O(t). 



(52) 



We will abuse notation slightly by using the same name 
for the function and its Fourier transform, the argument 
and context serving to distinguish them. Thus, R s (t) has 
the Fourier transform 



zeroes and poles only in the upper half-plane. We denote 
the roots of the numerator of P (as a function of ioS) as 



x± 



1 1 / 



(55) 



One possible choice of s is to make x+ > 0, X- < 0. Let 
us analyze this case first. This allows us to determine P± 
uniquely. 



iu) + 1 — s 

Now divide through Eq. (|54)l by P-(uj) to obtain 



(56) 



P+(u)h+(uj) = - 



ioj + 1 — s 



r_(w + is)(iu) + 1 — s) 



[iu — x_)(zu; — s) 



The first term on the right hand side has poles in both 
half-planes, and we separate them by partial fractions 



P+{u)h+{u) 

1 + X- — S 



(s — a;_ ) (iuj — s) 
. r_ (oj-\-is)(icu-\-l — s) 

(s - x-)(iu - x-) IUJ ~ X - ' 



(58) 



Since the product P+h + is guaranteed by construction 
to have poles only in the lower half-plane the terms with 
the poles in the upper half-plane on the right hand side 
must separately vanish, and we obtain 



(s — x—)(iu> — x+) 



(59) 



Going back to the fictitious "time" variable t and using 
g+(t) = e st h + (t), we obtain 



R s (u) = / dte luJt R s (t) = - 



1 — s —iu) - 



(53) 



The existence of the Fourier transform implies < s < 1 , 
but does not choose s uniquely. In general, the choice 
of s determines the class of functions which are allowed 
as solutions, as we will see explicitly below. Now the 
equation becomes 



(iu — s) (iu+l — s) 



h+(u) 



+ r_(w + is). 



(54) 



To proceed further we separate the prefactor of h+ on 
the left hand side [call it P(uj)] into a product P(u>) = 
P + (uj)P-(oj), where, by construction, P+(u) has zeroes 
and poles only in the lower half-plane and P-(u>) has 



<?+(*) 



0(t)e 



(s — X + )t 



(s-x-) 



(60) 



Note that with the definitions of x± , s drops out of this 
expression. Translating back to the original variables, we 

see that Foc(-|) 22 .It is easily verified that, 
as assumed in Eq. p5|) . the integral converges as k — > 0. 
Thus, the Wiener-Hopf method demonstrates the power 
law behavior of the solution to the integral equation. 

Another choice of s would be to make both x± > 0, 
yielding a solution which diverges more slowly as k — > 0. 
The general solution is a linear combination of both these 
solutions, and interestingly we see that the equation in 
its present form does not uniquely specify a particular 
combination. This ambiguity is lifted by introducing a 
lower cutoff fco in momentum, corresponding to a finite 
system size. We show below using an alternate method 
how this leads to a unique solution. 
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2. Solution of Equivalent Differential Equation 



Eq. (|43l) can also be solved with the model kernel by 
converting it into a differential equation. 
Differentiating Eq. (|4"3"]l we get 



F'(k) 



F(k') 



(61) 



Differentiating this equation, we get 



F " {k) J 2 [ h r dk '~ h (t) F(fc,) " TS t dk 'h F W) 



k 2 



F(k)h 



(62) 



but from Eq. (f6"Tj) , the first 2 terms in the square brackets 
are simply — | 2F ^ . Therefore the differential equation 
corresponding to Eq. (14*31 is 



F"{k) + ~F>(k) + ^F(k) = 0, 

which has general solutions of the form 
F(k) = A+k x+ +A_k x -, 



(63) 



(64) 



with k = k/A, A± = and 7 = VI - 2/3. The 

coefficients A± are determined by substituting Eq. [64] 
back into the integral equation. This results in power 
law behavior for k 3> fco, with exponent A+, which goes 
from real to complex when [3 exceeds 1/2. Moreover, 
M (q — > 0) may be evaluated, yielding 



M(0) = 



A 



vf1 + j- (3 + k$(-l + 7 + /3) 
This has poles for /3 > 1/2 when 



\/2f3 — 1 In fco = 2 arctan _^ — — h 27rn, 



(65) 



(66) 



with integer n and < arctan (x) < it. Note that the 
distance between poles vanishes logarithmically as fco 
0, as discussed above. Furthermore, for j3 > 1/2, Hq 
becomes ill-defined unless an infinitesimal imaginary part 
is introduced in /?, so that /3 = 1/2 becomes a branch 
point for M(0). We interpret this as the signal of a phase 
transition in the thermodynamic limit, since M(0) need 
not be real and positive beyond this point. 



V. DENSITY RESPONSE 

In this final section we return to another measurable 
quantity, the density response function. To be concrete 



we will use our result to compute the induced charge 
around a Coulomb impurity, which generates the poten- 
tial Ze/er, where e = |e|. Our procedure is to first solve 
Eq. (|35[) numerically, from which we compute the (static) 
irreducible polarizability 



Il(q) = % 



d 2 k ~ 



(2nf 



K/3 2 p iaa Ti3 1 i3 2 (k,q). 



(67) 



The density response function is then computed by 
an RPA surn^, except that instead of using the non- 
interacting polarizability we use our irreducible polariz- 
ability, which includes excitonic corrections via the ladder 
diagrams. The result of this takes the form 



D(q) = 



-U(q) 



l+n(q)Uc(qY 



(68) 



where Uc {q) = 2-k^vf /q is the Coulomb interaction. Fi- 
nally, the Fourier transform of the induced electron den- 
sity is given by 



Sn(q) = D(q)(f> eKt (q). 



(69) 



In the Coulomb impurity case, the external potential 
<?W(<?) =-ZU c {q)- 

To find numerical solutions to Eq. (|35[) . we discretize 

the allowed values of momentum k and q, and replace in- 
tegrations by sums over the grid of allowed momenta. In 
doing this, some subtleties arise. Since we cannot retain 
an infinite number of momentum points, we must confine 
the sums to a finite region, most conveniently taken to 
be square. If we use the MDF spectrum and wavefunc- 
tions (needed to construct K) in Eq. [35] in this "Brillouin 
zone" , the former will be periodic but not the latter. The 
discontinuity in wavefunctions leads to spurious oscilla- 
tions in the final result. In principle this can be overcome 
by simulating the system on a honeycomb lattice and us- 
ing the full tight-binding spectrum and wavefunctions for 
graphene. However, this is numerically costly and unnec- 
essary, because the low energy physics is almost entirely 
determined by the wavefunctions and spectra near the 
Dirac cones. Moveover, since Coulomb interactions are 
relatively weak at short wavelengths, one can neglect the 
intervalley scattering, so that it should be sufficient to 
consider only one Dirac cone, whereas a simulation of a 
honeycomb lattice would force us to include two due to 
fermion doubling 1 . 

As a compromise we consider models that have simpler 
bandstructures than graphene but still have a Dirac cone. 
One such model arises in the theory of the surface of a 
topological insulato r 25 i 26 , and has the form (Hvf = 1) 



Cn-\-x T" C n — C n _|_y -\- (l.C. 



(70) 



(71) 
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where 

h(k) = (sink x )a x + (shxk y )a v + (m + cosfc x + cosk y ) a z . 

This is a model defined on a square lattice, with each 
site supporting a two-component vector of localized or- 
bitals combined into annihilation operators c ra , and with 
the sum over n running through all the lattice sites. The 
crystal momentum k is measured in units of 1/a, with a 
the lattice constant. For m = 2, there is a single Dirac 
point at the center of the BZ, and the spinor structure 
in the vicinity of this point is the same as near the Dirac 
points in graphene. Thus we expect this model to repro- 
duce the low-energy behavior of graphene. We adopt this 
model for our numerical solution of Eq. Q35p . 

A second subtlety arises in the doped case. For exam- 
ple, when the chemical potential /x > 0, the quantity K 
in Eq. (|55|) takes the form 



9{e % - M ) - 0(e- - M )l <?+ 1& (%+ Ql (k + q] 



^ k k-\-q 

( £ k - ^)9hp2 (k)9a 2ai (k+_q) 

k k+q 
£ k + £ k + q 

where eg is the positive eigenvalue of h(k) and g a p[k) = 
(^cOafag )p> with rf~ being the eigenvectors of h(k) cor- 
responding to ±eg respectively. 

Because of the step functions, a naive numerical inte- 
gration by a discrete summation works poorly, because 
the function being integrated is not smooth on the scale 
of the grid. This problem may be overcome using the Tri- 
angular Linear Analytic (TLA) method 2 "^. We divide 
the square Brillouin zone into small squares, and each 
small square is further subdivided into two right trangles 
along one of the diagonals. Weights for the integrand can 
then be assigned at the corners of the triangles employ- 
ing the parameterization formulas in Ref. |28p^. Using 
this weighting scheme to approximate the integrals gives 
far better results than a naive lattice sum. 

Finally, we note that U(q) in Eq. (|55|) is better rep- 
resented by the RPA screened Coulomb interaction than 
the bare Coulomb interaction, i.e., 



U(q) = 



Uc(q) 



1 + n RPA (q)U C (q) 



(74) 



The irreducible RPA polarizability II RPA (q) has been cal- 
culated by a number of authors (see, e.g. Refs. Ho| and 
l3lt ) . For the case of undoped graphene there is no qual- 
itative difference between using screened or unscreened 
Coulomb interactions in the ladder rungs, because the 
functional forms of U(q) and Uc(q) are the same, and 



0.6 
0.5 



non-interacting^ 
ladder approx. 



N 

S 0.4 




0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 



FIG. 7: (Color online) Total induced electron number density 
divided by Z in the undoped case. 



only the effective value of j3 is renormalized. When /i =/= 0, 
however, the Fermi surface introduces a length scale into 
the problem, allowing genuine screening of the Coulomb 
interaction at long distances. This can be modeled by 
using a contact interaction on the ladder rung o 31 ' 32 al- 
though we find this introduces problems at large wavevec- 
tors, as we describe below. 

Fig. [7] shows the total induced electron number density 
in the undoped case, together with the RPA result 



,RPA 



(q = 0) IT RPA (g)£/c(<7) 



mi 

8 



z 



l + U^ A (q)Uc(q) 



and the non-interacting result 



Sn non (q = 0) 



n 



RPA 



71^ 



(75) 



(76) 



for comparison. It is clear that the non-interacting result 
exceeds 1 for (3 > 8/tt, while the RPA result approaches 1 
as {3 — v co. Results for the total induced charge using the 
RPA and the ladder approximation are not qualitatively 
different. The ladder approximation result is larger than 
the RPA result, and the difference increases with f3. 

For doped graphene, any charged impurity will induce 
an equal and opposite screening charge, both in the RPA 
and when exciton corrections are included. However, we 
find in the latter case a strong quantitative difference be- 
tween the two in the shape of the screening cloud. This 
is due to the effect of the exciton corrections on the ir- 
reducible polarizability in the doped case, illustrated in 
Fig.m For small /3, II is close to the RPA result for MDF 
as expected, except that the for q < 2kp there is a nega- 
tive slope, and for q close to the Brillouin zone boundary 
it is significantly below the RPA result for MDF, sim- 
ply because of the presence of the zone boundary. For 
larger /?, the curve is higher and the deviation from the 
RPA result for MDF is larger, and for the largest (3, an 
additional "hump" structure develops just below 2kp, as 
illustrated in the inset of Fig. [8] We will consider this 
extra structure in more detail below. 

The cusp in the density response function at q = 2kF is 
well-known to induce Friedel oscillations^ 3 -. The change 
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FIG. 8: (Color online) The continuous color curves are the 
irreducible polarizability 11(g) for the model Hamiltonian 
Eq. (170[l in the doped case, with RPA screened Coulomb in- 
teraction as the rungs of the ladders, for f) = 0.2, 0.4, . . . , 1.0, 
with higher curves corresponding to higher /3. The dashed 
black curve is the RPA result for MDF for comparison. The 
inset is a blow-up showing the developing structure just below 
q = 2k,F- For all the numerical results shown in this figure 
/i = 0.225hv F /a. 
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FIG. 9: (Color online) Induced electron number density 
(times r 2 /Z) as a function of distance r from the impurity 
in the doped case. On the far right side, from top to bottom, 
the curves correspond to /3 — 0.2, 0.4, . . . , 1.0. For compari- 
son, the dashed line is the RPA result with /3 = 0.2. 



in shape of the irreducible polarizability near 2Uf essen- 
tially deepens this cusp, leading to an enhancement of 
these oscillations. Fig. [9] shows the induced charge dis- 
tribution, illustrating this enhancement. We also note 
that with the exciton corrections included, the induced 
charge density falls off somewhat faster than the RPA 
result. 

We can also carry out the calculation with contact in- 
teractions as rungs of the ladders, i.e. U(q) = uq in 



FIG. 10: (Color online) Irreducible polarizability in the ladder 
approximation with contact interaction as the rungs of the 
ladders, for MDF in the doped case [Eq. ((78}]. uo = uokp/vF- 
The inset is the same plot with a larger range of q. 



Eq. ((55)) , which then becomes 



d 2 q' 



(77) 



K a 



2P27172 



■ 7172 



This is a considerable simplification relative to the equa- 
tion we had to solve for the RPA-screened Coulomb in- 
teraction; we now have a simple set of linear equations 
for Y a2 p 2 (q). The resulting irreducible polarizability is 



n(<?) 



i 



^n" PA ( g ) 



(78) 



This is plotted in Fig. [TU] for the doped case. We see 
that the simpler contact interaction does not enhance 
the cusp at q = 2k p, suggesting that the correct long- 
distance form of the screened Coulomb potential is an 
important ingredient in obtaining this behavior. Note 
also that because of the minus sign in the denominator 
of Eq. ((75)) and the monotonic increase of II RPA (<7) at 
large q, there is a pole at some sufficiently large q for any 
positive Mo, suggesting an instability in this model which 
is absent for the more realistic U(q). 

We can also obtain results for rungs with contact inter- 
actions numerically for our topological insulator surface 
model, Eq. ((70)) . illustrated in Fig. 11(a) Comparison 



((70)) . illustrated in Fig. |ll(a) 
between this and the nearly analytic results for MDF's 
allow us to assess which features may be introduced by 
going from the latter to the former. As we can see, the 
curves are very similar in overall scale to those when 
the interaction is the RPA screened Coulomb interac- 
tion (Fig. [5]). Note however that the deepening of the 
2kp cusp is absent in both the numerical result and the 
analytical one, suggesting that our numerical results are 
reasonably accurate at small q. However, we note that 
for the largest values of 0, extra structure near 2kF de- 
velops that appears analogous to what we found in the 
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FIG. 11: (Color online) Numerical result for the irreducible 
polarizability in the ladder approximation with contact inter- 
action as the rungs of the ladders, for the model Eq. (I70|l in 
the doped case. The insets are the blow-up around q = 2fcp. 
(a)tto = 0.2,0.4,0.6,0.8,0.9; (b)fio = 1.0. 



Coulomb case. This structure appears only in the result 
for the model Eq. (fTO)) . not in the analytical result for 
MDF's. It is thus reasonable to assume that the analo- 
gous structure in the RPA screened Coulomb interaction 
case is peculiar to Eq. (|70|) as well. It is interesting to 
speculate then that one may be able to distinguish the 
Dirac cone in the graphene system from that of at the 
surface of a topological insulator through such structure 
at very large j3. 

Finally, in Fig. 11(b) we illustrate that divergent 
behavior emerges when uq is sufficiently large, which 
evolves into a double pole from the "hump" structure 
below 2kp . This suggests the system becomes unsta- 
ble for contact interactions of sufficiently large magni- 
tude, a behavior that occurs as we observed above for 
any uq when MDF's are subject to contact interactions. 
(Note, however, in this case the instability sets in for 
q < 2kp, whereas for contact interactions the instability 
occurs at much larger q when the interaction is weak.) 
While this behavior is absent in the screened Coulomb 
interaction case, it is possible that at very small distances 



where the atomic orbital physics becomes relevant such a 
contact model becomes appropriate. Since this instabil- 
ity appears in the density-density response function this 
naively suggests that there is a phase transition into a 
charge-density wave. However, other transitions - spin 
or valley density waves, for example - may preempt this 
transition. We leave the nature of such an instability and 
its applicability to real graphene as open questions^. 



VI. CONCLUSION 

In this work we have investigated excitonic effects 
for graphene with Coulomb interactions, as modeled by 
massless Dirac fermions. We have shown that there is 
power law behavior in a general 4-leg vertex function in 
the particle-hole channel. The exponent becomes com- 
plex as the coupling constant j3 is increased above a crit- 
ical value. This is analogous to what happens in the 
problem of a single MDF interacting with a charged im- 
purity. This non-analytic behavior can be canceled away 
for certain combinations of the vertex function, and we 
find in particular that it is absent in the density-density 
response function. It is however retained in a sublattice 
antisymmetric response function. Although the power 
law behavior originates due to short length scale physics 
(close approaches of particle-hole pairs), it impacts the 
physics at large distances because of the absence of a 
length scale in the Hamiltonian. For finite size systems 
the transition appears to be one with broken chiral sym- 
metry, inducing a gap in the spectrum; however, this 
interpretation breaks down in the thermodynamic limit. 
We speculate that in this case the transition involves the 
formation of a mass gap which fluctuates among differ- 
ent possible forms, and is a precursor to a true broken 
symmetry state which emerges at still larger values of the 
coupling p. 

We have also calculated the density response in the lad- 
der approximation numerically using a simplified model 
Hamiltonian that occurs in the context of topological in- 
sulators, which has only one Dirac point and a square 
Brillouin zone. The calculation was carried out for both 
undoped and doped cases. In the latter case we com- 
pared results for RPA screened Coulomb interactions and 
contact interactions in the rungs of the ladders. While 
we expect the former interaction to be more realistic, 
both interactions in many respects give similar results. 
For Coulomb interactions, we find a strongly enhanced 
cusp in the irreducible polarizability at 2kp, which leads 
to much stronger Friedel oscillations than expected from 
the RPA. We also find a hump-like structure at stronger 
interaction scales just below 2k p. For contact interac- 
tions this hump evolves into poles with increasing inter- 
action strength, whereas no pole is seen for RPA screened 
Coulomb interactions in the range of j3 we have studied. 
We presented evidence that the extra structure is pe- 
culiar to our model Hamiltonian, suggesting it may be 
present at the surface of a topological insulator. Finally, 
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we note that with contact interactions, MDF's do con- 
tain a pole at larger q for any positive u , suggesting a 
short wavelength instability. 
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Appendix A 

In this Appendix, we discuss some details of the TLA 
method. This may be viewed as a two-dimensional ver- 
sion of the "tetrahedron method" 35 ' 36 which is widely 
applied for Brillouin zone integrations of three dimen- 
sional systems. As mentioned in Ref. [28l . when using 
linear interpolation, the surfaces e = s(x, y) (in this Ap- 
pendix x = k x and y = k y ) are straight lines and the 
parameterizations are easily given (with e being one of 
the parameters). There are two cases: e < £2 (but larger 
than £1) and e > £2 (but smaller than £3). The param- 
eterizations are given by Eqs. (14) and (16) in Ref. l28l . 
respectively. For convenience we reproduce them here. 



k =fc 3 + - — — (£3 - kx) 



£3 - £1 
e - £3 



(h - k 2 ) - 



_£3 - £2 
for e > £2, and 



6 — £1 -* -* 

k =fci H (fc 3 — fci) 



£3 ~ £1 



(Al) 



£3 - £1 
e — £1 



£2 - £1 



(*a - h) - 



£1 



£3 - £1 



(*3-fcl) 



(A2) 



for e < £2; in both cases ^ u ^ 1. 

Our goal is to express the integral over a basic triangle 
in terms of the values of £(x, y) and the integrand at the 
three corners of the triangle, i.e. 



f 3 
I(fi) = / d 2 k9(££ - n)f(k) « 22wi(£ 1 ,£ 2 ,£ 3 ,^)f(ki), 

(A3) 

where J A means integration over the triangle, Wi, i = 
1,2,3 are the weights at the three corners (remember 
that the corners are labeled so that £\ < £2 < £3). We 
will also use the shorthand /, = f(ki) below. 

There are four possibilities regarding the value of /i as 
compared to £i,£2,£3: 

(i) [i ^ £ 3 : 



(k 3 , £3) 



(k2, £2) 




(£3, £3) 



(k2, £2) 




(fci.ei) 



(a) 



(b) 



FIG. 12: (Color online) Illustrations for (a) case (ii) and 
(b) case (iii). In both subfigures the red and purple arrows 
label the directions of increasing e and u respectively. How- 
ever, while they correspond to Eq. (I All) in (a), in (b) they 
correspond to Eq. (IA2|) . 



This is trivial and the result isu^=0,i = l,2,3. 
(ii) £ 2 < M < £3: 

In this case, we need to use Eq. (|A1|) for the parame- 
terization, and 



/(/*)« / d 2 kf(k(e,u)) 



de / du 



2A(£ 3 - e) 
(£3 - £2) (£3 - £1) 



f(k(e,u)) 



(A4) 



where J means integration over the uppermost (light blue) triangle in Fig. 12(a) and the extra factor in the integrand 
is just the Jacobian, with A being the area of the triangle 123 (not the shaded triangle). 
We intepolate / linearly, i.e. 



f(x,y) Kipi +p 2 x+p 3 y = ^2pi9i(x,y), 



(A5) 



where g\ = 1, g 2 = x, g 3 = y, and the coefficients pi, i = 1,2,3 are determined by solving the equations 
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Y,i=iPi9i{xj,Vj) = fj- The result is 

{x 3 y 2 - x 2 y 3 )fi + (xiy 3 - x 3 y x )f 2 + (x 2 yi ~ x x y 2 )f 3 , » c s 

Pi = ; ; > (A6) 

x 2 yi - x 3 yi - xiy 2 + x 3 y 2 + x x y 3 - x 2 y 3 

p 2 and p 3 have the same denominator but with the numerator being [(y 3 — y 2 )fi + (y± — y 3 )f 2 + (y 2 — yi)f 3 ] for p 2 
and [{x 2 - x 3 )fi + (x 3 - x\)fa + (x x - x 2 )f 3 ] for p 3 . 

2(e 3 -e) 



g rS3 pi 

h {p) ~ J]] Api I de I du 

iu JO 



(£3 - e 2 )(e 3 - £1) 



The results of the integrations are 
2(£ 3 - e) 



(£3 - £2)^3 - ei) 
(£3 - £2)^3 - £1) 



du 



(£3 - m) 2 



(£3-£2)(£3-£i)' 



(A7) 



(A8) 



du{x 3 + — — -7- (2:3 — x\) + u 



(£3 - m) 2 



(£3 - £2)(£3 - £1) 



^3 - g(£3 - P) 



£3 - £l 

x 3 - xi x 3 - x 2 
£3 — £1 £3 ~ £2 



e - £3 , s e - e 3 

■(x 3 - x 2 ) - 



£3 - £2 



£3 - £1 



(23 - X\) 



(A9) 



and h(p) 3 is the same as Ii(/J.) 2 but with the cc's replaced by y's. Finally, we have 

A[(x 3 y 2 - x 2 y 3 )h{ij)i + (y 3 - y 2 )h{p) 2 + (x 2 - x 3 )Ii(n) a ] 

wi = ■ ■ , A10) 

x 2 yi - x 3 yi - xiy 2 + x 3 y 2 + x x y 3 - x 2 y 3 

w 2 and w 3 have the same denominator but with the numerator being A[(x\y 3 — x 3 yi)I\(n)\ + (y% — y 3 )h(p.) 2 + {x 3 — 
xi)h(p) 3 ] for w 2 and A[(x 2 yx - xiy 2 )h(p)i + (y 2 - yi)h{p) 2 + {x\ - x 2 )h(p,) 3 ] for w 3 . 
(iii)ei < n < £2 

In this case, I{p) ~ J d rkf(k(e,u)) where the integration is over the uppermost (cyan) triangle and the yellow 
quadrilateral in Fig. 12(b) The former is approximately the same (to linear order of the size of the basic triangles) 
as h(e 2 ), while the latter is 



3 ^ <.£ 2 2(g gj) ^ ^ 

Iu(n)&y~]Api de- ^— r / dug l = y^Apihi{p,) i . 

^ Jp (£2 ~£l)(£3 ~£l) Jo 



(All) 



The in the formula given above for case (ii) should be replaced by 



h{e 2 )i + hi{p)i 



h{s 2 ) 2 + ^11(^)2 



£3 ~ £2 
£3 - £1 



de 



2(e-ei) 



(£2 - £l)(£3 - El) Jo 



du 



£3 - £2 , (£2 - £i) 2 - (m~ £i) 2 



£3 -El (£2 ~ £l)(£3 - £1) 



(A12) 



£3 - £2 



£3 - £1 

du { x\ 



£3 - g(£3 - £2) 



X 







e — £1 
£3 - £1 



£3 ~ £2 
£3 - £1 



£3 - -(£3 - £2) 



x 3 - Xi x 3 - x 2 
£3 - £1 £3 — £2 
e — £1 
£2 - £1 
x 3 - xi x 3 - x 2 
£3 — £1 £3 — £2 
1 



de- 



2(e 



(x 3 -Xi)+U 



(x 2 - Xl) 



(e 2 - £i)(£ 3 - £1) 
e — £1 



£3 - £1 
1 



-(x 3 - Xi) 



(£2 - £l)(£3 - £1) 



[(£2 - £i) 2 - (M - £i) 2 ] 2:1 + 3 [(£2 ~ £i) 3 - (M - £i) 3 ] 



X2 — X1 x 3 — xi 



£2 — £1 £3 — £1 

and iii(£2) 3 + J[i(m)3) which is the same as ^(£2)2 + hi(n) 2 but with the x's replaced by y's. 
(xv)fi < £1 

This case is trivial too, we simply have w x — w 2 — w 3 = 1/3. 
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